rm(list = ls())
## llamando a las librerias que (creemos que) vamos a usar
options("install.lock"=FALSE)
if (!require ("pacman")) install.packages ("pacman")

pacman::p_load (dplyr
                , haven
                , ggplot2
                , ggpubr
                , ggthemes  ## si quieren mas themes
                , readstata13
                , readxl
                , sf
                , tidyverse
                , tidyr
                , units
                , viridis ## paleta de colores Viridis
                , wesanderson## p/usar paleta de colores de Wes Anderson
                , stringr
                , RColorBrewer
                , patchwork
                , Rmisc
                , lfe
                , stargazer
                , AER
                , haven
                , skimr
                , modelsummary
                , terra
                , fixest
                , vtable
                , did
                , cowplot
                , grid)

setwd("/Users/mateovasquez/Dropbox/")

#####

base <- read.dta13("master_violence_landtitle_caracteristicas.dta")

rm(list = ls())
## llamando a las librerias que (creemos que) vamos a usar
options("install.lock"=FALSE)
if (!require ("pacman")) install.packages ("pacman")

pacman::p_load (dplyr
                , haven
                , ggplot2
                , ggpubr
                , ggthemes  ## si quieren mas themes
                , readstata13
                , readxl
                , sf
                , tidyverse
                , tidyr
                , units
                , viridis ## paleta de colores Viridis
                , wesanderson## p/usar paleta de colores de Wes Anderson
                , stringr
                , RColorBrewer
                , patchwork
                , Rmisc
                , lfe
                , stargazer
                , AER
                , haven
                , skimr
                , modelsummary
                , terra
                , fixest
                , vtable
                , did
                , cowplot
                , grid)

setwd("/Users/mateovasquez/Dropbox/")

#####

base <- read.dta13("master_violence_landtitle_caracteristicas.dta")

#Outcome Variable
names(base)

base <- base %>% 
  dplyr::mutate(total_violence = rowSums(across(c(sumattacks_public, sumattacks_guerrilla,sumattacks_paramilitary)), na.rm = TRUE),
                health_inst = rowSums(across(c(health_post, health_center,hospital)), na.rm = TRUE))


#Renaming Variables 

base <- base %>% 
  dplyr::rename("indian_title_area_cs" = "area_resguardos_indigenas_cs",
                "indian_title_area" = "area_resguardos_indigenas",
                "afro_title_area_cs" = "area_titul_comunidades_negras_cs",
                "afro_title_area" = "area_titul_comunidades_negras")

#Treatment 	

base <- base %>% 
  dplyr::mutate(post = ifelse(Año>1995,1,0),
                any_titleXpost = any_title*post,
                any_title2 = any_title,
                post2 = post)


#the mean level of violence in neighboring municipalities (1988–1995)
Neighbors <- base %>% 
  dplyr::filter(Año<=1995) %>% 
  dplyr::group_by(codmpio) %>% 
  dplyr::summarise(spatlagguerrsum100_100_1995 = mean(spatlagguerrsum100_100, na.rm = T))

summary(Neighbors$spatlagguerrsum100_100_1995)

base <- left_join(base,Neighbors,
                  by = "codmpio")

#Controls X post
base <- base %>%
  dplyr::mutate(
    altura_post = altura*post, #ELEVATION
    rainfall_post = rainfall*post, #RAINFALL
    coca_0_1_post = coca_0_1*post, # COCA
    presence_mines1560_post = presence_mines1560*post, #GOLD MINES
    num_encomiendas1560_post = num_encomiendas1560*post, #ENCOMIENDAS
    royalroaddist_post = royalroaddist*post, #distance to royal roads
    distancia_mercado_post = distancia_mercado*post, #nearest market towns (km)
    pop95_post = pop95*post, #municipal population in 1995 
    pct_minority85_post = pct_minority85*post, #share of minority population in 1985 (%)
    conflictos_1901_1917_post = conflictos_1901_1917*post,#historical conflict using data on the inci- dence of land conflict between 1901 and 1931
    conflictos_1918_1931_post = conflictos_1918_1931*post,#historical conflict using data on the inci- dence of land conflict between 1901 and 1931
    Violencia_48_a_53_post = Violencia_48_a_53*post, #presence of armed combat during La Violencia (1948–1953)
    anucraids1971_78_post = anucraids1971_78*post, #frequency of land invasions between 1971 and 1978
    spatlagguerrsum100_100_1995_post = spatlagguerrsum100_100_1995*post, #the mean level of violence in neighboring municipalities (1988–1995)
    ltotal_plots_rfmd_1960_85_post = ltotal_plots_rfmd_1960_85*post, # total number of plots reformed between 1960 and 1985 to proxy for rural grievances.
    totviolencia85_post = totviolencia85*post) #ojo: La Violencia (1948 - 1958)

# Applying asinh transformation to each specified variable directly in mutate

base <- base %>%
  mutate(
    as_total_violence = asinh(total_violence),
    as_sumattacks_public = asinh(sumattacks_public),
    as_sumattacks_paramilitary = asinh(sumattacks_paramilitary),
    as_sumattacks_guerrilla = asinh(sumattacks_guerrilla))

# Running the regression with felm, including fixed effects and clustering

base <- base %>%
  mutate(Depto_year = paste(coddepto, Año, sep = "_"))


####PARAMILITARY GRAPH
twfe.est <- feols(as_sumattacks_paramilitary ~ i(Año,factor(any_title2), 1995) | 
                    codmpio,  
                  data = base, cluster = "codmpio")

twfe.est.output <- cbind.data.frame(twfe.est$coefficients,
                                    SE = twfe.est$se) %>% 
  rownames_to_column()

twfe.est.output <- as.data.frame(twfe.est.output)

twfe.est.output$Year <- substr(twfe.est.output$rowname, 6, 9)
twfe.est.output$Group <- substr(twfe.est.output$rowname, 31, 31)

twfe.est.output <- twfe.est.output %>% 
  dplyr::rename("Estimate" = "twfe.est$coefficients") %>% 
  dplyr::mutate(Year = as.numeric(Year),
                LB_CI = Estimate-(1.96*SE),
                UB_CI = Estimate+(1.96*SE),
                Model = "TWFE",
                interval = ifelse(Year < 1985, 1980,
                                  ifelse(Year < 1989,1985,
                                         ifelse(Year < 1994,1990,
                                                ifelse(Year < 2000,1995,
                                                       ifelse(Year < 2005,2000,
                                                              ifelse(Year < 2010,2005,2010))))))) %>% 
  dplyr::filter(!is.na(Year)) %>% 
  dplyr::select(-c(SE,rowname))

table(twfe.est.output$Year, twfe.est.output$interval)

results <- twfe.est.output %>%
  dplyr::group_by(interval, Group) %>%
  dplyr::summarise(Estimate_mean = mean(Estimate, na.rm = TRUE),
                   LB_CI_mean = mean(LB_CI, na.rm = TRUE),
                   UB_CI_mean = mean(UB_CI, na.rm = TRUE))

# Assuming 'results' is your data frame and 'Group' is the column you're adjusting
results$Group <- factor(results$Group, levels = c("1", "0"))

# Then, your plotting code follows
Paramilitary <- ggplot(results, aes(x = interval, y = Estimate_mean, shape = as.factor(Group))) +
  geom_line(aes(linetype=as.factor(Group)), size = 1) +
  scale_linetype_manual(values=c("solid","dashed"),
                        breaks=c("1", "0"), # Ensure this matches the new order
                        labels=c("Titled (Law 70)", "Not Titled")) + # Update labels to match the new order
  geom_point(size = 3) +
  scale_shape_manual(values = c("1" = 18, "0" = 15),
                     breaks=c("1", "0"),
                     labels=c("Titled (Law 70)", "Not Titled")) +
  geom_errorbar(aes(ymin=LB_CI_mean, ymax=UB_CI_mean, group= factor(Group)), size=0.2, width=0.2) +
  labs(x = "Interval", y = "Marginal Effect Coefficient") +
  geom_hline(yintercept = mean(base$as_sumattacks_paramilitary), col = "red") +
  geom_vline(xintercept = 1995, col = "black") +
  theme_calc() +
  theme(legend.position = c(0.1, 0.9),
        legend.justification = c(0, 1),
        legend.text = element_text(size = 12),
        legend.key.size = unit(1.5, "lines"),
        legend.background = element_rect(fill = "white", colour = "black")) +
  scale_x_continuous(breaks = seq(1980, 2010, by = 5), limits = c(1980, 2010)) +
  scale_y_continuous(breaks = seq(-.4, .6, by = .1), limits = c(-.5, .7)) +
  ggtitle("Paramilitary Attacks") +
  guides(shape = guide_legend(title = NULL), linetype = guide_legend(title = NULL))

Paramilitary

#####POLICE/ARMY ATTACKS

twfe.est <- feols(as_sumattacks_public ~ i(Año,factor(any_title2), 1995) | 
                    codmpio,  
                  data = base, cluster = "codmpio")

twfe.est.output <- cbind.data.frame(twfe.est$coefficients,
                                    SE = twfe.est$se) %>% 
  rownames_to_column()

twfe.est.output <- as.data.frame(twfe.est.output)

twfe.est.output$Year <- substr(twfe.est.output$rowname, 6, 9)
twfe.est.output$Group <- substr(twfe.est.output$rowname, 31, 31)

twfe.est.output <- twfe.est.output %>% 
  dplyr::rename("Estimate" = "twfe.est$coefficients") %>% 
  dplyr::mutate(Year = as.numeric(Year),
                LB_CI = Estimate-(1.96*SE),
                UB_CI = Estimate+(1.96*SE),
                Model = "TWFE",
                interval = ifelse(Year < 1985, 1980,
                                  ifelse(Year < 1989,1985,
                                         ifelse(Year < 1994,1990,
                                                ifelse(Year < 2000,1995,
                                                       ifelse(Year < 2005,2000,
                                                              ifelse(Year < 2010,2005,2010))))))) %>% 
  dplyr::filter(!is.na(Year)) %>% 
  dplyr::select(-c(SE,rowname))

table(twfe.est.output$Year, twfe.est.output$interval)

results <- twfe.est.output %>%
  dplyr::group_by(interval, Group) %>%
  dplyr::summarise(Estimate_mean = mean(Estimate, na.rm = TRUE),
                   LB_CI_mean = mean(LB_CI, na.rm = TRUE),
                   UB_CI_mean = mean(UB_CI, na.rm = TRUE))

# Assuming 'results' is your data frame and 'Group' is the column you're adjusting
results$Group <- factor(results$Group, levels = c("1", "0"))

# Then, your plotting code follows
public <- ggplot(results, aes(x = interval, y = Estimate_mean, shape = as.factor(Group))) +
  geom_line(aes(linetype=as.factor(Group)), size = 1) +
  scale_linetype_manual(values=c("solid","dashed"),
                        breaks=c("1", "0"), # Ensure this matches the new order
                        labels=c("Titled (Law 70)", "Not Titled")) + # Update labels to match the new order
  geom_point(size = 3) +
  scale_shape_manual(values = c("1" = 18, "0" = 15),
                     breaks=c("1", "0"),
                     labels=c("Titled (Law 70)", "Not Titled")) +
  geom_errorbar(aes(ymin=LB_CI_mean, ymax=UB_CI_mean, group= factor(Group)), size=0.2, width=0.2) +
  labs(x = "Interval", y = "Marginal Effect Coefficient") +
  geom_hline(yintercept = mean(base$as_sumattacks_public), col = "red") +
  geom_vline(xintercept = 1995, col = "black") +
  theme_calc() +
  theme(legend.position = "none") +
  scale_x_continuous(breaks = seq(1980, 2010, by = 5), limits = c(1980, 2010)) +
  scale_y_continuous(breaks = seq(-.4, .6, by = .1), limits = c(-.5, .63)) +
  ggtitle("Public Attacks") +
  guides(shape = guide_legend(title = NULL), linetype = guide_legend(title = NULL))

public

####GUERRILLA

twfe.est <- feols(as_sumattacks_guerrilla ~ i(Año,factor(any_title2), 1995) | 
                    codmpio,  
                  data = base, cluster = "codmpio")

twfe.est.output <- cbind.data.frame(twfe.est$coefficients,
                                    SE = twfe.est$se) %>% 
  rownames_to_column()

twfe.est.output <- as.data.frame(twfe.est.output)

twfe.est.output$Year <- substr(twfe.est.output$rowname, 6, 9)
twfe.est.output$Group <- substr(twfe.est.output$rowname, 31, 31)

twfe.est.output <- twfe.est.output %>% 
  dplyr::rename("Estimate" = "twfe.est$coefficients") %>% 
  dplyr::mutate(Year = as.numeric(Year),
                LB_CI = Estimate-(1.96*SE),
                UB_CI = Estimate+(1.96*SE),
                Model = "TWFE",
                interval = ifelse(Year < 1985, 1980,
                                  ifelse(Year < 1989,1985,
                                         ifelse(Year < 1994,1990,
                                                ifelse(Year < 2000,1995,
                                                       ifelse(Year < 2005,2000,
                                                              ifelse(Year < 2010,2005,2010))))))) %>% 
  dplyr::filter(!is.na(Year)) %>% 
  dplyr::select(-c(SE,rowname))

table(twfe.est.output$Year, twfe.est.output$interval)

results <- twfe.est.output %>%
  dplyr::group_by(interval, Group) %>%
  dplyr::summarise(Estimate_mean = mean(Estimate, na.rm = TRUE),
                   LB_CI_mean = mean(LB_CI, na.rm = TRUE),
                   UB_CI_mean = mean(UB_CI, na.rm = TRUE))

# Assuming 'results' is your data frame and 'Group' is the column you're adjusting
results$Group <- factor(results$Group, levels = c("1", "0"))

# Then, your plotting code follows
guerrilla <- ggplot(results, aes(x = interval, y = Estimate_mean, shape = as.factor(Group))) +
  geom_line(aes(linetype=as.factor(Group)), size = 1) +
  scale_linetype_manual(values=c("solid","dashed"),
                        breaks=c("1", "0"), # Ensure this matches the new order
                        labels=c("Titled (Law 70)", "Not Titled")) + # Update labels to match the new order
  geom_point(size = 3) +
  scale_shape_manual(values = c("1" = 18, "0" = 15),
                     breaks=c("1", "0"),
                     labels=c("Titled (Law 70)", "Not Titled")) +
  geom_errorbar(aes(ymin=LB_CI_mean, ymax=UB_CI_mean, group= factor(Group)), size=0.2, width=0.2) +
  labs(x = "Interval", y = "Marginal Effect Coefficient") +
  geom_hline(yintercept = mean(base$as_sumattacks_guerrilla), col = "red") +
  geom_vline(xintercept = 1995, col = "black") +
  theme_calc() +
  theme(legend.position = "none") +
  scale_x_continuous(breaks = seq(1980, 2010, by = 5), limits = c(1980, 2010)) +
  scale_y_continuous(breaks = seq(-.4, .6, by = .1), limits = c(-.5, .63)) +
  ggtitle("Guerrilla Attacks") +
  guides(shape = guide_legend(title = NULL), linetype = guide_legend(title = NULL))

guerrilla

Paramilitary + public + guerrilla

#####
#pacific
####PARAMILITARY GRAPH
twfe.est <- feols(as_sumattacks_paramilitary ~ i(Año,factor(any_title2), 1995) | 
                    codmpio,  
                  data = base %>% dplyr::filter(gpacifica == 1), cluster = "codmpio")

twfe.est.output <- cbind.data.frame(twfe.est$coefficients,
                                    SE = twfe.est$se) %>% 
  rownames_to_column()

twfe.est.output <- as.data.frame(twfe.est.output)

twfe.est.output$Year <- substr(twfe.est.output$rowname, 6, 9)
twfe.est.output$Group <- substr(twfe.est.output$rowname, 31, 31)

twfe.est.output <- twfe.est.output %>% 
  dplyr::rename("Estimate" = "twfe.est$coefficients") %>% 
  dplyr::mutate(Year = as.numeric(Year),
                LB_CI = Estimate-(1.96*SE),
                UB_CI = Estimate+(1.96*SE),
                Model = "TWFE",
                interval = ifelse(Year < 1985, 1980,
                                  ifelse(Year < 1989,1985,
                                         ifelse(Year < 1994,1990,
                                                ifelse(Year < 2000,1995,
                                                       ifelse(Year < 2005,2000,
                                                              ifelse(Year < 2010,2005,2010))))))) %>% 
  dplyr::filter(!is.na(Year)) %>% 
  dplyr::select(-c(SE,rowname))

table(twfe.est.output$Year, twfe.est.output$interval)

results <- twfe.est.output %>%
  dplyr::group_by(interval, Group) %>%
  dplyr::summarise(Estimate_mean = mean(Estimate, na.rm = TRUE),
                   LB_CI_mean = mean(LB_CI, na.rm = TRUE),
                   UB_CI_mean = mean(UB_CI, na.rm = TRUE))

# Assuming 'results' is your data frame and 'Group' is the column you're adjusting
results$Group <- factor(results$Group, levels = c("1", "0"))

# Then, your plotting code follows
Paramilitary_pac <- ggplot(results, aes(x = interval, y = Estimate_mean, shape = as.factor(Group))) +
  geom_line(aes(linetype=as.factor(Group)), size = 1) +
  scale_linetype_manual(values=c("solid","dashed"),
                        breaks=c("1", "0"), # Ensure this matches the new order
                        labels=c("Titled (Law 70)", "Not Titled")) + # Update labels to match the new order
  geom_point(size = 3) +
  scale_shape_manual(values = c("1" = 18, "0" = 15),
                     breaks=c("1", "0"),
                     labels=c("Titled (Law 70)", "Not Titled")) +
  geom_errorbar(aes(ymin=LB_CI_mean, ymax=UB_CI_mean, group= factor(Group)), size=0.2, width=0.2) +
  labs(x = "Interval", y = "Marginal Effect Coefficient") +
  geom_hline(yintercept = mean(base$as_sumattacks_paramilitary[base$gpacifica==1]), col = "red") +
  geom_vline(xintercept = 1995, col = "black") +
  theme_calc() +
  theme(legend.position = c(0.1, 0.9),
        legend.justification = c(0, 1),
        legend.text = element_text(size = 12),
        legend.key.size = unit(1.5, "lines"),
        legend.background = element_rect(fill = "white", colour = "black")) +
  scale_x_continuous(breaks = seq(1980, 2010, by = 5), limits = c(1980, 2010)) +
  scale_y_continuous(breaks = seq(-.2, 1, by = .2), limits = c(-.2, 1.2)) +
  ggtitle("Paramilitary Attacks") +
  guides(shape = guide_legend(title = NULL), linetype = guide_legend(title = NULL))

Paramilitary_pac

#####POLICE/ARMY ATTACKS

twfe.est <- feols(as_sumattacks_public ~ i(Año,factor(any_title2), 1995) | 
                    codmpio,  
                  data = base %>% dplyr::filter(gpacifica == 1), cluster = "codmpio")

twfe.est.output <- cbind.data.frame(twfe.est$coefficients,
                                    SE = twfe.est$se) %>% 
  rownames_to_column()

twfe.est.output <- as.data.frame(twfe.est.output)

twfe.est.output$Year <- substr(twfe.est.output$rowname, 6, 9)
twfe.est.output$Group <- substr(twfe.est.output$rowname, 31, 31)

twfe.est.output <- twfe.est.output %>% 
  dplyr::rename("Estimate" = "twfe.est$coefficients") %>% 
  dplyr::mutate(Year = as.numeric(Year),
                LB_CI = Estimate-(1.96*SE),
                UB_CI = Estimate+(1.96*SE),
                Model = "TWFE",
                interval = ifelse(Year < 1985, 1980,
                                  ifelse(Year < 1989,1985,
                                         ifelse(Year < 1994,1990,
                                                ifelse(Year < 2000,1995,
                                                       ifelse(Year < 2005,2000,
                                                              ifelse(Year < 2010,2005,2010))))))) %>% 
  dplyr::filter(!is.na(Year)) %>% 
  dplyr::select(-c(SE,rowname))

table(twfe.est.output$Year, twfe.est.output$interval)

results <- twfe.est.output %>%
  dplyr::group_by(interval, Group) %>%
  dplyr::summarise(Estimate_mean = mean(Estimate, na.rm = TRUE),
                   LB_CI_mean = mean(LB_CI, na.rm = TRUE),
                   UB_CI_mean = mean(UB_CI, na.rm = TRUE))

# Assuming 'results' is your data frame and 'Group' is the column you're adjusting
results$Group <- factor(results$Group, levels = c("1", "0"))

# Then, your plotting code follows
public_pac <- ggplot(results, aes(x = interval, y = Estimate_mean, shape = as.factor(Group))) +
  geom_line(aes(linetype=as.factor(Group)), size = 1) +
  scale_linetype_manual(values=c("solid","dashed"),
                        breaks=c("1", "0"), # Ensure this matches the new order
                        labels=c("Titled (Law 70)", "Not Titled")) + # Update labels to match the new order
  geom_point(size = 3) +
  scale_shape_manual(values = c("1" = 18, "0" = 15),
                     breaks=c("1", "0"),
                     labels=c("Titled (Law 70)", "Not Titled")) +
  geom_errorbar(aes(ymin=LB_CI_mean, ymax=UB_CI_mean, group= factor(Group)), size=0.2, width=0.2) +
  labs(x = "Interval", y = "Marginal Effect Coefficient") +
  geom_hline(yintercept = mean(base$as_sumattacks_public[base$gpacifica==1]), col = "red") +
  geom_vline(xintercept = 1995, col = "black") +
  theme_calc() +
  theme(legend.position = "none") +
  scale_x_continuous(breaks = seq(1980, 2010, by = 5), limits = c(1980, 2010)) +
  scale_y_continuous(breaks = seq(-.2, 1, by = .2), limits = c(-.2, 1.2)) +
  ggtitle("Public Attacks") +
  guides(shape = guide_legend(title = NULL), linetype = guide_legend(title = NULL))

public_pac

####GUERRILLA

twfe.est <- feols(as_sumattacks_guerrilla ~ i(Año,factor(any_title2), 1995) | 
                    codmpio,  
                  data = base %>% dplyr::filter(gpacifica == 1), cluster = "codmpio")

twfe.est.output <- cbind.data.frame(twfe.est$coefficients,
                                    SE = twfe.est$se) %>% 
  rownames_to_column()

twfe.est.output <- as.data.frame(twfe.est.output)

twfe.est.output$Year <- substr(twfe.est.output$rowname, 6, 9)
twfe.est.output$Group <- substr(twfe.est.output$rowname, 31, 31)

twfe.est.output <- twfe.est.output %>% 
  dplyr::rename("Estimate" = "twfe.est$coefficients") %>% 
  dplyr::mutate(Year = as.numeric(Year),
                LB_CI = Estimate-(1.96*SE),
                UB_CI = Estimate+(1.96*SE),
                Model = "TWFE",
                interval = ifelse(Year < 1985, 1980,
                                  ifelse(Year < 1989,1985,
                                         ifelse(Year < 1994,1990,
                                                ifelse(Year < 2000,1995,
                                                       ifelse(Year < 2005,2000,
                                                              ifelse(Year < 2010,2005,2010))))))) %>% 
  dplyr::filter(!is.na(Year)) %>% 
  dplyr::select(-c(SE,rowname))

table(twfe.est.output$Year, twfe.est.output$interval)

results <- twfe.est.output %>%
  dplyr::group_by(interval, Group) %>%
  dplyr::summarise(Estimate_mean = mean(Estimate, na.rm = TRUE),
                   LB_CI_mean = mean(LB_CI, na.rm = TRUE),
                   UB_CI_mean = mean(UB_CI, na.rm = TRUE))

# Assuming 'results' is your data frame and 'Group' is the column you're adjusting
results$Group <- factor(results$Group, levels = c("1", "0"))

# Then, your plotting code follows
guerrilla_pac <- ggplot(results, aes(x = interval, y = Estimate_mean, shape = as.factor(Group))) +
  geom_line(aes(linetype=as.factor(Group)), size = 1) +
  scale_linetype_manual(values=c("solid","dashed"),
                        breaks=c("1", "0"), # Ensure this matches the new order
                        labels=c("Titled (Law 70)", "Not Titled")) + # Update labels to match the new order
  geom_point(size = 3) +
  scale_shape_manual(values = c("1" = 18, "0" = 15),
                     breaks=c("1", "0"),
                     labels=c("Titled (Law 70)", "Not Titled")) +
  geom_errorbar(aes(ymin=LB_CI_mean, ymax=UB_CI_mean, group= factor(Group)), size=0.2, width=0.2) +
  labs(x = "Interval", y = "Marginal Effect Coefficient") +
  geom_hline(yintercept = mean(base$as_sumattacks_guerrilla[base$gpacifica==1]), col = "red") +
  geom_vline(xintercept = 1995, col = "black") +
  theme_calc() +
  theme(legend.position = "none") +
  scale_x_continuous(breaks = seq(1980, 2010, by = 5), limits = c(1980, 2010)) +
  scale_y_continuous(breaks = seq(-.2, 1, by = .2), limits = c(-.2, 1.2)) +
  ggtitle("Guerrilla Attacks") +
  guides(shape = guide_legend(title = NULL), linetype = guide_legend(title = NULL))

guerrilla_pac

todas <- (Paramilitary + public + guerrilla)/(Paramilitary_pac + public_pac + guerrilla_pac)

